import intake
from easygems import healpix as egh
import healpy as hp

import matplotlib.pyplot as plt
import warnings

warnings.filterwarnings("ignore", category=FutureWarning)

current_location = "NCAR" #"online"
cat = intake.open_catalog("https://digital-earths-global-hackathon.github.io/catalog/catalog.yaml")[current_location]
#list(cat)

scream = cat['scream2D_hrly'](zoom = 7).to_dask()
scream_z8 = cat['scream2D_hrly'](zoom = 8).to_dask()
scream_z9 = cat['scream2D_hrly'](zoom = 9).to_dask()
scream_z10 = cat['scream2D_hrly'](zoom = 10).to_dask()

ERA5

rda_url           = 'https://data.rda.ucar.edu/'
era5_catalog      = rda_url + 'pythia_era5_24/pythia_intake_catalogs/era5_catalog.json'
#rda_data          = "/gpfs/csfs1/collections/rda/data/"
#era5_catalog      = rda_data + 'pythia_era5_24/pythia_intake_catalogs/era5_catalog.json'
era5_cat = intake.open_esm_datastore(era5_catalog)
#era5_cat.df[['long_name','variable', 'units']].drop_duplicates()[190:]#.head()
pre_mean_cat = era5_cat.search(variable='MTPR',frequency = 'hourly',year=2020,month=6)
pre_cat = era5_cat.search(variable='MXTPR',frequency = 'hourly',year=2020,month=6)
#pre_cat.df['path']

/glade/u/apps/opt/conda/envs/2025-digital-earths-global-hackathon/lib/python3.12/site-packages/intake_esm/cat.py:249: DtypeWarning: Columns (2) have mixed types. Specify dtype option on import or set low_memory=False.
  df = pd.read_csv(

xarray_open_kwargs = {
    'engine': 'h5netcdf',
    'chunks': {},  # Specify any chunking if needed
    'backend_kwargs': {}  # Any additional backend arguments if required
}
ds_pre_era5 = pre_cat.to_dataset_dict(xarray_open_kwargs=xarray_open_kwargs)
ds_pre_era5.keys()
ds_pre_mean_era5 = pre_mean_cat.to_dataset_dict(xarray_open_kwargs=xarray_open_kwargs)

--> The keys in the returned dictionary of datasets are constructed as follows:
	'datatype.step_type'

100.00% [1/1 00:03<00:00]

--> The keys in the returned dictionary of datasets are constructed as follows:
	'datatype.step_type'

100.00% [1/1 00:04<00:00]

pre_mean_era5 = ds_pre_mean_era5['fc.meanflux'].MTPR*3.6e3
pre_era5 = ds_pre_era5['fc.minmax'].MXTPR
pre_era5
pre_era5_mmhr = pre_era5*3.6e3
lc = pre_era5_mmhr.coords["longitude"]
la = pre_era5_mmhr.coords["latitude"]
pre_era5_june_subset = pre_era5_mmhr.loc[dict(longitude=lc[(lc > 105) & (lc < 122)], latitude=la[(la > 25) & (la < 33)])]

pre_era5_june=pre_era5_june_subset.mean(axis=(0,1))
pre_era5_june

<xarray.DataArray 'MXTPR' (latitude: 31, longitude: 67)> Size: 8kB
dask.array<mean_agg-aggregate, shape=(31, 67), dtype=float32, chunksize=(31, 67), chunktype=numpy.ndarray>
Coordinates:
  * latitude   (latitude) float64 248B 32.75 32.5 32.25 ... 25.75 25.5 25.25
  * longitude  (longitude) float64 536B 105.2 105.5 105.8 ... 121.2 121.5 121.8

xarray.DataArray

'MXTPR'

• latitude: 31
• longitude: 67

• dask.array<chunksize=(31, 67), meta=np.ndarray>

  Array Chunk Bytes 8.11 kiB 8.11 kiB Shape (31, 67) (31, 67) Dask graph 1 chunks in 15 graph layers Data type float32 numpy.ndarray

• Coordinates: (2)
  • latitude
    (latitude)
    float64
    32.75 32.5 32.25 ... 25.5 25.25

    long_name :

    latitude

    short_name :

    lat

    units :

    degrees_north

    array([32.75, 32.5 , 32.25, 32.  , 31.75, 31.5 , 31.25, 31.  , 30.75, 30.5 ,
           30.25, 30.  , 29.75, 29.5 , 29.25, 29.  , 28.75, 28.5 , 28.25, 28.  ,
           27.75, 27.5 , 27.25, 27.  , 26.75, 26.5 , 26.25, 26.  , 25.75, 25.5 ,
           25.25])

  • longitude
    (longitude)
    float64
    105.2 105.5 105.8 ... 121.5 121.8

    long_name :

    longitude

    short_name :

    lon

    units :

    degrees_east

    array([105.25, 105.5 , 105.75, 106.  , 106.25, 106.5 , 106.75, 107.  , 107.25,
           107.5 , 107.75, 108.  , 108.25, 108.5 , 108.75, 109.  , 109.25, 109.5 ,
           109.75, 110.  , 110.25, 110.5 , 110.75, 111.  , 111.25, 111.5 , 111.75,
           112.  , 112.25, 112.5 , 112.75, 113.  , 113.25, 113.5 , 113.75, 114.  ,
           114.25, 114.5 , 114.75, 115.  , 115.25, 115.5 , 115.75, 116.  , 116.25,
           116.5 , 116.75, 117.  , 117.25, 117.5 , 117.75, 118.  , 118.25, 118.5 ,
           118.75, 119.  , 119.25, 119.5 , 119.75, 120.  , 120.25, 120.5 , 120.75,
           121.  , 121.25, 121.5 , 121.75])

• Indexes: (2)
  • latitude
    PandasIndex

    PandasIndex(Index([32.75,  32.5, 32.25,  32.0, 31.75,  31.5, 31.25,  31.0, 30.75,  30.5,
           30.25,  30.0, 29.75,  29.5, 29.25,  29.0, 28.75,  28.5, 28.25,  28.0,
           27.75,  27.5, 27.25,  27.0, 26.75,  26.5, 26.25,  26.0, 25.75,  25.5,
           25.25],
          dtype='float64', name='latitude'))

  • longitude
    PandasIndex

    PandasIndex(Index([105.25,  105.5, 105.75,  106.0, 106.25,  106.5, 106.75,  107.0, 107.25,
            107.5, 107.75,  108.0, 108.25,  108.5, 108.75,  109.0, 109.25,  109.5,
           109.75,  110.0, 110.25,  110.5, 110.75,  111.0, 111.25,  111.5, 111.75,
            112.0, 112.25,  112.5, 112.75,  113.0, 113.25,  113.5, 113.75,  114.0,
           114.25,  114.5, 114.75,  115.0, 115.25,  115.5, 115.75,  116.0, 116.25,
            116.5, 116.75,  117.0, 117.25,  117.5, 117.75,  118.0, 118.25,  118.5,
           118.75,  119.0, 119.25,  119.5, 119.75,  120.0, 120.25,  120.5, 120.75,
            121.0, 121.25,  121.5, 121.75],
          dtype='float64', name='longitude'))

• Attributes: (0)

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cf
import cmocean

projection = ccrs.Robinson(central_longitude=180)
fig, axis = plt.subplots(1, 1, subplot_kw=dict(projection=projection))
#pre_era5_june.plot(
pre_mean_era5.mean(axis=(0,1)).plot(
    ax=axis,
    cmap=cmocean.cm.thermal,
    transform=ccrs.PlateCarree(),
    levels=[0.1,0.2,0.3,0.4,0.5,0.6]
)
extent = [90, 124, 20, 55]
axis.set_extent(extent, crs=ccrs.PlateCarree())
axis.coastlines()
axis.add_feature(cf.BORDERS)
axis.set_title("ERA5 Precip (mm/hr)", fontsize = 18 )

Text(0.5, 1.0, 'ERA5 Precip (mm/hr)')

SCREAM

import cftime 
start_time = cftime.DatetimeNoLeap(2020, 6, 1, 0, 0, 0, 0)
end_time = cftime.DatetimeNoLeap(2020, 6, 30, 0, 0, 0, 0)

import uxarray as ux

uxds = ux.UxDataset.from_healpix(scream)
scream_precip = uxds.pr
#scream_precip
scream_precip_mmhr = scream_precip * 36e5
scream_precip_mmhr_june =  scream_precip_mmhr.sel(time=slice(start_time, end_time))

uxds = ux.UxDataset.from_healpix(scream_z8)
scream8_precip = uxds.pr
scream8_precip_mmhr = scream8_precip * 36e5
scream8_precip_mmhr_june =  scream8_precip_mmhr.sel(time=slice(start_time, end_time))

uxds = ux.UxDataset.from_healpix(scream_z9)
scream9_precip = uxds.pr
scream9_precip_mmhr = scream9_precip * 36e5
scream9_precip_mmhr_june =  scream9_precip_mmhr.sel(time=slice(start_time, end_time))

uxds = ux.UxDataset.from_healpix(scream_z10)
scream10_precip = uxds.pr
scream10_precip_mmhr = scream10_precip * 36e5
scream10_precip_mmhr_june =  scream10_precip_mmhr.sel(time=slice(start_time, end_time))

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cf
import cmocean
from easygems import healpix as egh

projection = ccrs.Robinson(central_longitude=180)
'''
scream_precip_mmhr_june.mean(dim="time").compute().plot(
    projection=projection,
    cmap=cmocean.cm.thermal,
    features=["borders", "coastline"],
    title="SCREAM Precip (zoom=7) June 2020",
    width=1100,
    clim=(0.1,0.6)
)
'''
fig, axis = plt.subplots(1, 1, subplot_kw=dict(projection=projection))
extent = [90, 124, 20, 55]
axis.set_extent(extent, crs=ccrs.PlateCarree())
#egh.healpix_show(scream_precip_june, ax=axes[0], cmap=cmp, vmin=.1, vmax=.6)
axis.add_feature(cf.COASTLINE, linewidth=0.8)
axis.add_feature(cf.BORDERS, linewidth=0.4)
axis.set_title("SCREAM Precip (mm/hr)", fontsize = 18 )
fg = egh.healpix_show(scream_precip_mmhr_june.mean(dim="time"), ax=axis, cmap=cmocean.cm.thermal, vmin=.1, vmax=.6)#cmocean.cm.thermal, vmax=0.5)
fig.colorbar(fg)
plt.show()

lon_bounds = (105, 122)
lat_bounds = (25, 33)

scream_precip_june_subset = scream_precip_mmhr_june.subset.bounding_box(lon_bounds, lat_bounds)

scream8_precip_june_subset = scream8_precip_mmhr_june.subset.bounding_box(lon_bounds, lat_bounds)

scream9_precip_june_subset = scream9_precip_mmhr_june.subset.bounding_box(lon_bounds, lat_bounds)

scream10_precip_june_subset = scream10_precip_mmhr_june.subset.bounding_box(lon_bounds, lat_bounds)

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cf
import cmocean

projection = ccrs.Robinson(central_longitude=180)

scream_precip_june_subset.mean(dim="time").compute().plot(
    projection=projection,
    cmap=cmocean.cm.thermal,
    features=["borders", "coastline"],
    title="SCREAM Precip (zoom=7) June 2020",
    width=1100,
    clim=(0.1,0.6)
)

WARNING:param.GeoOverlayPlot01111: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.

import numpy as np
bins = np.logspace(np.log10(0.5), np.log10(30), 300)
histogram_scream, _ = np.histogram(
    scream_precip_june_subset, bins, density=True
)
histogram_era5, _ = np.histogram(
    pre_era5_mmhr_subset, bins, density=True)

histogram_scream8, _ = np.histogram(
    scream8_precip_june_subset, bins, density=True)
histogram_scream9, _ = np.histogram(
    scream9_precip_june_subset, bins, density=True)
histogram_scream10, _ = np.histogram(
    scream10_precip_june_subset, bins, density=True)

import seaborn as sns
xaxis = (bins[1:] + bins[:-1]) / 2
fig, ax = plt.subplots(figsize=(8, 6), facecolor="white")
plt.plot(xaxis, histogram_scream, label="SCREAM (zoom=7)")#"Yangtze River Region, CHN")
plt.plot(xaxis, histogram_era5, label="ERA5")
plt.plot(xaxis, histogram_scream8, label="SCREAM (zoom=8)")
plt.plot(xaxis, histogram_scream9, label="SCREAM (zoom=9)")
plt.plot(xaxis, histogram_scream10, label="SCREAM (zoom=10)")
ax.set_yscale("log")
ax.set_xlabel("P / mmhr$^{-1}$")
ax.set_ylabel("PDF")
sns.despine(offset=10)
plt.legend()

<matplotlib.legend.Legend at 0x146e76c45820>